/*
 * fastaAlg.cpp
 *
 *  Created on: 30-09-2011
 *      Author: pawel
 */

#include"fastaAlg.h"
#include"../algorithms.h"
#include<pthread.h>
#include"../smart2DArray.h"
#include"../seqFeederList.h"
#include"../scoreMatrix.h"
#include"fastaGraph.h"
#include"../stat/cStat.h"
#include<Bpp/Seq/Sequence.h>
#include<string.h>
#include<sys/param.h>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/dijkstra_shortest_paths.hpp>
#include <boost/graph/dag_shortest_paths.hpp>

#include<assert.h>
#include<stdio.h>
#include<iostream>

fastaAlg::fastaAlg(const optAlgThreadParam &optAlgParam,size_t heurThreads,size_t optThreads,const fastaAlgParam &fastaParam)
:heurAlgV(optAlgParam,heurThreads,optThreads),fastaParam(fastaParam),qLength((optAlgThParam.querySeqPtr->getContent()).size()),
smartQueryString(1,qLength),smartQueryContent(1,qLength)
{
	char *queryCString = (smartQueryString.getPtr())[0];
	int *contentCTable = (smartQueryContent.getPtr())[0];
	memcpy(queryCString,(optAlgThParam.querySeqPtr->toString().c_str()),qLength * sizeof(char));
	memcpy(contentCTable,&(optAlgParam.querySeqPtr->getContent()[0]),qLength * sizeof(int));
	smart2DArray<char> smartKTupleString(1,fastaParam.k_tuple + 1);
	char *kTupleCString = (smartKTupleString.getPtr())[0];
	kTupleCString[fastaParam.k_tuple] = 0;
	for(size_t i = 0;i<=qLength - fastaParam.k_tuple;i++){
		memcpy(kTupleCString,(queryCString + i),fastaParam.k_tuple * sizeof(char) );
		queryKTouples.insert(std::pair<std::string,int>(kTupleCString,i));

	}



}
void fastaAlg::runAlg(){
	heurAlgV::startThreads(this);
}
void* fastaAlg::runThread(){
	bpp::BasicSequence *dbSequencePtr;
	size_t lic=0;
	while(true){
		//Consumer
		RESTART:
		pthread_mutex_lock( mainSeqFeederPtr->getMutex());
			if(mainSeqFeederPtr->isClosed()){
				pthread_mutex_unlock( mainSeqFeederPtr->getMutex());
				break;
			}
		while(mainSeqFeederPtr->isEmpty()){
			if(mainSeqFeederPtr->isProductionEnded()){
				mainSeqFeederPtr->setClose();
				pthread_mutex_unlock( mainSeqFeederPtr->getMutex() );
				goto RESTART;
			}
			pthread_cond_wait(mainSeqFeederPtr->getNotEmpty(), mainSeqFeederPtr->getMutex());
		}
		try{
			dbSequencePtr = mainSeqFeederPtr->getSeq();


		}catch(std::exception& e){
			pthread_mutex_unlock( mainSeqFeederPtr->getMutex() );
			pthread_cond_signal (mainSeqFeederPtr->getNotFull());
			continue;
		}
		pthread_mutex_unlock( mainSeqFeederPtr->getMutex() );
		pthread_cond_signal (mainSeqFeederPtr->getNotFull());
		bounds optAlgBounds;
		{
			__gnu_cxx::slist<diagonal> diags;
			diags = computeDiagonalHits(dbSequencePtr);
			rescoreDiagonals(diags,dbSequencePtr);
			optAlgBounds = computeBounds(diags);


		}
		//Producer
		pthread_mutex_lock(optAlgThParam.seqFeederPtr->getMutex());
		while(optAlgThParam.seqFeederPtr->isFull()){
			pthread_cond_wait(optAlgThParam.seqFeederPtr->getNotFull(),optAlgThParam.seqFeederPtr->getMutex());
		}
		optAlgThParam.seqFeederPtr->pushSeq(dbSequencePtr);
		seqAlignParams addParams(&optAlgBounds);
		optAlgThParam.seqFeederPtr->pushAddParams(&addParams,dbSequencePtr->getName());
		pthread_mutex_unlock(optAlgThParam.seqFeederPtr->getMutex());
		pthread_cond_signal (optAlgThParam.seqFeederPtr->getNotEmpty());

	}
	return 0;
}
__gnu_cxx::slist<diagonal> fastaAlg::computeDiagonalHits(bpp::BasicSequence *dbSeqPtr){

	__gnu_cxx::slist<diagonal> bestDiagonals;
	size_t seqLen = dbSeqPtr->size();
	smart2DArray<char> dbSmartSeqString(1,seqLen);
	char *dbSeqCString = (dbSmartSeqString.getPtr())[0];
	memcpy(dbSeqCString,dbSeqPtr->toString().c_str(),seqLen * sizeof(char));
	std::pair<boost::unordered_multimap<std::string,int>::iterator,boost::unordered_multimap<std::string,int>::iterator> matches;
	boost::unordered_multimap<std::string,int>::iterator iter;
	smart2DArray<char> smartKTupleStr(1,fastaParam.k_tuple + 1);
	char *kTupleStr = (smartKTupleStr.getPtr())[0];
	kTupleStr[fastaParam.k_tuple] = 0;
	smart2DArray<diagonal> smartDiagonalTab(1,qLength + seqLen);
	diagonal initDiag;
	initDiag.count=0;
	initDiag.diagNum=-qLength;
	initDiag.bestFragment.begin=0;
	initDiag.bestFragment.end=0;
	initDiag.bestFragment.score=0;
	diagonal *diagonalCTab = (smartDiagonalTab.getPtr())[0];
	for(size_t i=0;i<qLength + seqLen ;i++){
		diagonalCTab[i]= initDiag;
		initDiag.diagNum++;
	}
	size_t hitDiags=0;
	if(seqLen > fastaParam.k_tuple){
		for(size_t i = 0;i<= seqLen - fastaParam.k_tuple;i++)
		{
			memcpy(kTupleStr,(dbSeqCString + i),fastaParam.k_tuple * sizeof(char));
			matches = queryKTouples.equal_range(kTupleStr);
			for(iter =  matches.first; iter != matches.second;++iter){
				(diagonalCTab[ i -(*iter).second    + qLength]).count++;

			}
		}



		qsort(diagonalCTab,qLength + seqLen ,sizeof(diagonal),diagonal::cntCompare);
		statResult stats;
		stats = computeStats(diagonalCTab,qLength + seqLen);
		size_t tresh = round(stats.mean + 2.0* stats.stdDev );
		size_t bestDiags =10;
		bestDiags = MIN(bestDiags,qLength + seqLen );

		for(size_t i=0;i<bestDiags;i++){
			if(diagonalCTab[i].count == 0)break;

			bestDiagonals.push_front(diagonalCTab[i]);
			hitDiags ++;
		}

	}

	return bestDiagonals;
}
int fastaAlg::rescoreDiagonals(__gnu_cxx::slist<diagonal> &list,bpp::BasicSequence *dbSeqPtr){
	size_t dbContLen = dbSeqPtr->size();
	smart2DArray<int> smartdbContent(1,dbContLen);
	int *dbContentCArray = (smartdbContent.getPtr())[0];
	memcpy(dbContentCArray,&(dbSeqPtr->getContent()[0]),dbContLen * sizeof(int));
	char **scoreMatrix  = optAlgThParam.scoreMatrixPtr->getMatrix();
	int *queryContent = (smartQueryContent.getPtr())[0];
	__gnu_cxx::slist<diagonal>::iterator iter;
	__gnu_cxx::slist<diagonal>ndiags;
	int diagNum=0;
	size_t dbContCnt=0;
	size_t queryContCnt=0;
	for(iter = list.begin(); iter != list.end();++iter){
		diagNum = (*iter).diagNum;
		dbContCnt =0;
		queryContCnt=0;
		if(diagNum >0){
			dbContCnt = diagNum;
		}
		else if(diagNum < 0){
			queryContCnt = -diagNum;
		}
		int score;
		int scoreSum=0;
		(*iter).bestFragment.score = 0;
		int cnt=0;
		int cnt2=0;
		int best= INT_MIN;
		while(queryContCnt < qLength && dbContCnt < dbContLen){
			score = scoreMatrix[queryContent[queryContCnt]][dbContentCArray[dbContCnt]];
			scoreSum+= score;
			if(scoreSum > (*iter).bestFragment.score ){
				(*iter).bestFragment.score = scoreSum;
				(*iter).bestFragment.begin = cnt;
				(*iter).bestFragment.end = cnt2;
			}
			else  {
				if((*iter).bestFragment.begin -(*iter).bestFragment.end >2 && (*iter).bestFragment.score > 1 ){

					ndiags.push_front((*iter));
				}
				(*iter).bestFragment.score = 0;
				cnt = cnt2 + 1;
				scoreSum =0;
			}

			queryContCnt++;
			dbContCnt++;
			cnt2++;
		}

	}


	for(iter = ndiags.begin(); iter != ndiags.end();++iter){
		list.push_front(*iter);
	}


	return 0;
}
using namespace boost;
bounds fastaAlg::computeBounds(__gnu_cxx::slist<diagonal> &list){
	 typedef adjacency_list<vecS, vecS, directedS,
	    property<vertex_distance_t, int>, property<edge_weight_t, int> > graph_t;
	__gnu_cxx::slist<diagonal>::iterator iter1,iter2;
	size_t nodes = list.size();
	graph_t alignGraph(nodes);
	property_map<graph_t, vertex_distance_t>::type d_map = get(vertex_distance, alignGraph);
	smart2DArray<diagonal> smartDiagonalTab(1,nodes);
	diagonal *diagonalTab = (smartDiagonalTab.getPtr())[0];
	smart2DArray<size_t> smartNextVertex(1,nodes);
	size_t count=0;

	for(iter1 = list.begin();iter1 != list.end();++iter1){
		diagonalTab[count]=*(iter1);

		count++;
	}
	if(count == 0){
		bounds noHitsBounds;
		noHitsBounds.maxDiag = INT_MAX;
		noHitsBounds.minDiag = INT_MIN;
		return noHitsBounds;

	}

	bool edgesAdded=false;
	matrixPoint endI,beginJ,minEnd;
	size_t gapLength=0;
	int weight=0;
	size_t beginPointPosition=0;

	minEnd.i = UINT_MAX;
	minEnd.j = UINT_MAX;

	for(size_t i=0;i<nodes;i++){
		endI = diagonal::diagToEndPoint(diagonalTab[i]);
		size_t minLen = UINT_MAX;
		for(size_t j=0;j<nodes;j++){
			if(i == j)continue;

			beginJ = diagonal::diagToBegPoint(diagonalTab[j]);



			if(endI.i < beginJ.i && endI.j < beginJ.j ){
				gapLength = abs(diagonalTab[i].diagNum - diagonalTab[j].diagNum );
				{
					size_t distance = (abs(endI.i - beginJ.i) + abs(endI.j - beginJ.j));
					weight = -(diagonalTab[i].bestFragment.score + diagonalTab[j].bestFragment.score -( abs(endI.i - beginJ.i) + abs(endI.j - beginJ.j)) );

					{
						add_edge(i, j,weight , alignGraph);
						edgesAdded = true;

					}

				}
			}

		}
	}
	bounds retBounds;
	retBounds.maxDiag = INT_MIN;
	retBounds.minDiag = INT_MAX;
	if(edgesAdded){
		int minScoreG = INT_MAX;
		do{
			property_map<graph_t, vertex_distance_t>::type d_map = get(vertex_distance, alignGraph);
			std::vector<default_color_type> color(num_vertices(alignGraph));
			  std::vector<std::size_t> pred(num_vertices(alignGraph));
			  default_dijkstra_visitor vis;
			  std::less<int> compare;
			  closed_plus<int> combine;
			  property_map<graph_t, edge_weight_t>::type w_map = get(edge_weight, alignGraph);
			  dag_shortest_paths(alignGraph, beginPointPosition, d_map, w_map, &color[0], &pred[0],vis, compare, combine, (std::numeric_limits<int>::max)(), 0);
			  graph_traits<graph_t>::vertex_iterator vi , vi_end;
			  int minScore=INT_MAX;
			  size_t minLengthVertex=0;
			   for (boost::tie(vi, vi_end) = vertices(alignGraph); vi != vi_end; ++vi)
			   {
				   if(d_map[*vi] <= minScore && *vi != beginPointPosition){
						 minScore = d_map[*vi];
						minLengthVertex = *vi;
				   }
			   }

			   size_t currentVertex;
			   size_t nextVertex = minLengthVertex;

			   if(minScore < minScoreG){
				   retBounds.maxDiag = INT_MIN;
				   retBounds.minDiag = INT_MAX;
				   minScoreG = minScore;

				   do{

					   currentVertex = nextVertex;
						if(diagonalTab[currentVertex].diagNum > retBounds.maxDiag){
							retBounds.maxDiag = diagonalTab[currentVertex].diagNum;
						}
						if(diagonalTab[currentVertex].diagNum < retBounds.minDiag){
							retBounds.minDiag = diagonalTab[currentVertex].diagNum;
						}
					   nextVertex = pred[currentVertex];


				   }while(nextVertex != currentVertex);


			   }
			   beginPointPosition++;

		}while(beginPointPosition < nodes);
	}
	else
	{
		int bscore = INT_MIN;
		int bdiag=0;
		for(size_t i =0;i<nodes;i++){
			if(diagonalTab[i].bestFragment.score > bscore){
				bscore = diagonalTab[i].bestFragment.score;
				bdiag = diagonalTab[i].diagNum;
			}
		}

		retBounds.maxDiag = bdiag;
		retBounds.minDiag = bdiag;

	}
	if(retBounds.minDiag != INT_MAX){retBounds.minDiag-=fastaParam.w;}
	else{retBounds.minDiag = INT_MIN;}
	if(retBounds.maxDiag != INT_MIN){retBounds.maxDiag+=fastaParam.w;}
	else{retBounds.maxDiag = INT_MAX;}


return retBounds;
}
fastaAlg::~fastaAlg(){

}


